Contents

CESM2-Large Ensemble Reproduction of Kay et al. 2015

Introduction

This Jupyter Notebook demonstrates how one might use the NCAR Community Earth System Model v2 (CESM2) Large Ensemble (LENS) data hosted on AWS S3. The notebook shows how to reproduce figures 2 and 4 from the Kay et al. (2015) paper describing the CESM LENS dataset (doi:10.1175/BAMS-D-13-00255.1), with the LENS2 dataset.

This resource is intended to be helpful for people not familiar with elements of the Pangeo framework including Jupyter Notebooks, Xarray, and Zarr data format, or with the original paper, so it includes additional explanation.

Imports

%matplotlib inline
import warnings

warnings.filterwarnings("ignore")
# Silence dask.distributed logs
import logging

logger = logging.getLogger("distributed.utils_perf")
logger.setLevel(logging.ERROR)

import intake
import numpy as np
import pandas as pd
import xarray as xr
import hvplot.pandas, hvplot.xarray
import holoviews as hv
from ncar_jobqueue import NCARCluster
from distributed import Client
hv.extension('bokeh')

Spin up a Cluster

cluster = NCARCluster(memory="50GB", walltime='2:00:00', cores=4, processes=1, resource_spec='select=1:ncpus=1:mem=50GB')
cluster.scale(20)
client = Client(cluster)
client

Client

Client-7310d807-2633-11ec-848e-3cecef1b11fa

Connection method: Cluster object Cluster type: dask_jobqueue.PBSCluster
Dashboard: https://jupyterhub.hpc.ucar.edu/stable/user/mgrover/proxy/37785/status

Cluster Info

Read in Data

Use the Intake-ESM Catalog to Access the Data

catalog = intake.open_esm_datastore(
    '../data/glade-cesm2-le.json'
)
catalog

aws-cesm2-le catalog with 23 dataset(s) from 205 asset(s):

unique
component 3
spatial_domain 2
frequency 3
path 196
experiment 2
forcing_variant 2
variable 35
long_name 35
start_time 4
end_time 4

Subset for Daily Temperature Data (TREFHT)

We use Intake-ESM here to query for our dataset, focusing the smoothed biomass burning (smbb) experiment!

catalog_subset = catalog.search(variable='TREFHT', frequency='daily')
catalog_subset

aws-cesm2-le catalog with 4 dataset(s) from 4 asset(s):

unique
component 1
spatial_domain 1
frequency 1
path 4
experiment 2
forcing_variant 2
variable 1
long_name 1
start_time 2
end_time 2 
catalog_subset.df
component spatial_domain frequency path experiment forcing_variant variable long_name start_time end_time
0 atm global daily /glade/scratch/mgrover/data/lens2-aws/atm/dail... historical cmip6 TREFHT Reference height temperature 1850-01-01 12:00:00 2015-01-01 12:00:00
1 atm global daily /glade/scratch/mgrover/data/lens2-aws/atm/dail... historical smbb TREFHT Reference height temperature 1850-01-01 12:00:00 2015-01-01 12:00:00
2 atm global daily /glade/scratch/mgrover/data/lens2-aws/atm/dail... ssp370 cmip6 TREFHT Reference height temperature 2015-01-01 12:00:00 2101-01-01 12:00:00
3 atm global daily /glade/scratch/mgrover/data/lens2-aws/atm/dail... ssp370 smbb TREFHT Reference height temperature 2015-01-01 12:00:00 2101-01-01 12:00:00

Taking a look at the dataframe, we see there are two files - one for the historical run, and the other a future scenario (ssp370)

catalog_subset.df
component spatial_domain frequency path experiment forcing_variant variable long_name start_time end_time
0 atm global daily /glade/scratch/mgrover/data/lens2-aws/atm/dail... historical cmip6 TREFHT Reference height temperature 1850-01-01 12:00:00 2015-01-01 12:00:00
1 atm global daily /glade/scratch/mgrover/data/lens2-aws/atm/dail... historical smbb TREFHT Reference height temperature 1850-01-01 12:00:00 2015-01-01 12:00:00
2 atm global daily /glade/scratch/mgrover/data/lens2-aws/atm/dail... ssp370 cmip6 TREFHT Reference height temperature 2015-01-01 12:00:00 2101-01-01 12:00:00
3 atm global daily /glade/scratch/mgrover/data/lens2-aws/atm/dail... ssp370 smbb TREFHT Reference height temperature 2015-01-01 12:00:00 2101-01-01 12:00:00

Load in our datasets using .to_dataset_dict()

We use .to_dataset_dict() to return a dictionary of datasets, which we can now use for the analysis

dsets = catalog_subset.to_dataset_dict()

--> The keys in the returned dictionary of datasets are constructed as follows:
	'component.experiment.frequency.forcing_variant'
100.00% [4/4 00:07<00:00]

Here are the keys for each dataset - you will notice they follow the groupby attributes, with component.experiment.frequency.forcing_variant, one for the future and the other historical here

dsets.keys()

dict_keys(['atm.ssp370.daily.smbb', 'atm.historical.daily.smbb', 'atm.historical.daily.cmip6', 'atm.ssp370.daily.cmip6'])

Let’s load these into separate datasets, then merge these together!

historical_smbb = dsets['atm.historical.daily.smbb']
future_smbb = dsets['atm.ssp370.daily.smbb']

historical_cmip6 = dsets['atm.historical.daily.cmip6']
future_cmip6 = dsets['atm.ssp370.daily.cmip6']

merge_ds_smbb = xr.concat([historical_smbb, future_smbb], dim='time')
merge_ds_smbb = merge_ds_smbb.dropna(dim='member_id')

merge_ds_cmip6= xr.concat([historical_cmip6, future_cmip6], dim='time')
merge_ds_cmip6 = merge_ds_cmip6.dropna(dim='member_id')

merge_ds_smbb['time'] = merge_ds_smbb.indexes['time'].to_datetimeindex()
merge_ds_cmip6['time'] = merge_ds_cmip6.indexes['time'].to_datetimeindex()

Finally, we can load in just the variable we are interested in, Reference height temperature (TREFHT). We will separate this into two arrays:

  • t_smbb - the entire period of data for the smoothed biomass burning experiment (1850 to 2100)

  • t_cmip6 - the entire period of data for the cmip6 biomass burning experiment (1850 to 2100)

  • t_ref - the reference period used in the Kay et al. 2015 paper (1961 to 1990)

t_smbb = merge_ds_smbb.TREFHT
t_cmip6 = merge_ds_cmip6.TREFHT
t_ref = t_cmip6.sel(time=slice('1961', '1990'))

Read in the Grid Data

We also have Zarr stores of grid data, such as the area of each grid cell. We will need to follow a similar process here, querying for our experiment and extracting the dataset

grid_subset = catalog.search(component='atm', frequency='static', experiment='historical', forcing_variant='smbb')

_, grid = grid_subset.to_dataset_dict(aggregate=False, zarr_kwargs={"consolidated": True}).popitem()
grid

--> The keys in the returned dictionary of datasets are constructed as follows:
	'component.spatial_domain.frequency.path.experiment.forcing_variant.variable.long_name.start_time.end_time'
100.00% [1/1 00:00<00:00] 
<xarray.Dataset>
Dimensions:   (lat: 192, lon: 288, ilev: 31, lev: 30, bnds: 2, slat: 191, slon: 288)
Coordinates: (12/20)
    P0        float64 ...
    area      (lat, lon) float32 dask.array<chunksize=(192, 288), meta=np.ndarray>
    gw        (lat) float64 dask.array<chunksize=(192,), meta=np.ndarray>
    hyai      (ilev) float64 dask.array<chunksize=(31,), meta=np.ndarray>
    hyam      (lev) float64 dask.array<chunksize=(30,), meta=np.ndarray>
    hybi      (ilev) float64 dask.array<chunksize=(31,), meta=np.ndarray>
    ...        ...
    ntrm      int32 ...
    ntrn      int32 ...
  * slat      (slat) float64 -89.53 -88.59 -87.64 -86.7 ... 87.64 88.59 89.53
  * slon      (slon) float64 -0.625 0.625 1.875 3.125 ... 355.6 356.9 358.1
    w_stag    (slat) float64 dask.array<chunksize=(191,), meta=np.ndarray>
    wnummax   (lat) int32 dask.array<chunksize=(192,), meta=np.ndarray>
Dimensions without coordinates: bnds
Data variables:
    *empty*
Attributes:
    intake_esm_varname:      None
    intake_esm_dataset_key:  atm.global.static./glade/scratch/mgrover/data/le...

Extract the Grid Area and Compute the Total Area

Here, we extract the grid area and compute the total area across all grid cells - this will help when computing the weights…

ds_grid = xr.open_dataset('/glade/campaign/cgd/cesm/CESM2-LE/timeseries/atm/proc/tseries/month_1/AREA/b.e21.BHISTcmip6.f09_g17.LE2-1001.001.cam.h0.AREA.185001-185912.nc')

cell_area = ds_grid.isel(time=0).AREA.load()
total_area = cell_area.sum()
cell_area

<xarray.DataArray 'AREA' (lat: 192, lon: 288)>
array([[2.9948368e+07, 2.9948368e+07, 2.9948368e+07, ..., 2.9948368e+07,
        2.9948368e+07, 2.9948368e+07],
       [2.3957478e+08, 2.3957478e+08, 2.3957478e+08, ..., 2.3957478e+08,
        2.3957478e+08, 2.3957478e+08],
       [4.7908477e+08, 4.7908477e+08, 4.7908477e+08, ..., 4.7908477e+08,
        4.7908477e+08, 4.7908477e+08],
       ...,
       [4.7908477e+08, 4.7908477e+08, 4.7908477e+08, ..., 4.7908477e+08,
        4.7908477e+08, 4.7908477e+08],
       [2.3957478e+08, 2.3957478e+08, 2.3957478e+08, ..., 2.3957478e+08,
        2.3957478e+08, 2.3957478e+08],
       [2.9948368e+07, 2.9948368e+07, 2.9948368e+07, ..., 2.9948368e+07,
        2.9948368e+07, 2.9948368e+07]], dtype=float32)
Coordinates:
  * lat      (lat) float64 -90.0 -89.06 -88.12 -87.17 ... 87.17 88.12 89.06 90.0
  * lon      (lon) float64 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8
    time     object 1850-02-01 00:00:00
Attributes:
    units:         m2
    long_name:     area of grid box
    cell_methods:  time: mean

Compute the Weighted Annual Averages

Here, we utilize the resample(time="AS") method which does an annual resampling based on start of calendar year.

Note we also weight by the cell area, multiplying each value by the corresponding cell area, summing over the grid (lat, lon), then dividing by the total area

Setup the Computation - Below, we lazily prepare the calculation!

You’ll notice how quickly this cell runs - which is suspicious 👀 - we didn’t actually do any computation, but rather prepared the calculation

%%time
t_ref_ts = (
    (t_ref.resample(time="AS").mean("time")).weighted(cell_area).mean(('lat', 'lon'))
).mean(dim=("time", "member_id"))

t_smbb_ts = (
    (t_smbb.resample(time="AS").mean("time")).weighted(cell_area).mean(('lat', 'lon')))

t_cmip6_ts = (
    (t_cmip6.resample(time="AS").mean("time")).weighted(cell_area).mean(('lat', 'lon')))

CPU times: user 1.95 s, sys: 14.7 ms, total: 1.96 s
Wall time: 1.95 s

t_ref_ts.values

array(287.47504, dtype=float32)

Compute the Averages

At this point, we actually compute the averages - using %%time to help us measure how long it takes to compute…

%%time
t_ref_mean = t_ref_ts.compute()

CPU times: user 15.3 s, sys: 514 ms, total: 15.8 s
Wall time: 20.6 s

%%time
t_smbb_mean = t_smbb_ts.compute()

CPU times: user 1min 42s, sys: 3.94 s, total: 1min 46s
Wall time: 2min 3s

%%time
t_cmip6_mean = t_cmip6_ts.compute()

CPU times: user 2min 11s, sys: 4.99 s, total: 2min 16s
Wall time: 2min 42s

Convert to Dataframes

Our output data format is still an xarray.DataArray which isn’t neccessarily ideal… we could convert this to a dataframe which would be easier to work with!

We can convert our xarray.DataArray to a pandas dataframe using the following:

t_smbb_ts_df = t_smbb_mean.to_series().unstack().T
t_cmip6_ts_df = t_cmip6_mean.to_series().unstack().T

anomaly_smbb = (t_smbb_ts_df - t_ref_mean.data)
anomaly_cmip6 = (t_cmip6_ts_df - t_ref_mean.data)

Grab some Observational Data

In this case, we use the HADCRUT4 dataset

# Observational time series data for comparison with ensemble average
obsDataURL = "https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/cru/hadcrut4/air.mon.anom.median.nc"

ds = xr.open_dataset(obsDataURL).load()
ds

<xarray.Dataset>
Dimensions:    (lat: 36, lon: 72, time: 2060, nbnds: 2)
Coordinates:
  * lat        (lat) float32 87.5 82.5 77.5 72.5 ... -72.5 -77.5 -82.5 -87.5
  * lon        (lon) float32 -177.5 -172.5 -167.5 -162.5 ... 167.5 172.5 177.5
  * time       (time) datetime64[ns] 1850-01-01 1850-02-01 ... 2021-08-01
Dimensions without coordinates: nbnds
Data variables:
    time_bnds  (time, nbnds) datetime64[ns] 1850-01-01 1850-01-31 ... 2021-08-31
    air        (time, lat, lon) float32 nan nan nan nan nan ... nan nan nan nan
Attributes:
    platform:                        Surface
    title:                           HADCRUT4 Combined Air Temperature/SST An...
    history:                         Originally created at NOAA/ESRL PSD by C...
    Conventions:                     CF-1.0
    Comment:                         This dataset supersedes V3
    Source:                          Obtained from http://hadobs.metoffice.co...
    version:                         4.2.0
    dataset_title:                   HadCRUT4
    References:                      https://www.psl.noaa.gov/data/gridded/da...
    DODS_EXTRA.Unlimited_Dimension:  time

Compute the Weighted Temporal Mean from Seasons

def weighted_temporal_mean(ds):
    """
    weight by days in each month
    """
    time_bound_diff = ds.time_bnds.diff(dim="nbnds")[:, 0]
    wgts = time_bound_diff.groupby("time.year") / time_bound_diff.groupby(
        "time.year"
    ).sum(xr.ALL_DIMS)
    np.testing.assert_allclose(wgts.groupby("time.year").sum(xr.ALL_DIMS), 1.0)
    obs = ds["air"]
    cond = obs.isnull()
    ones = xr.where(cond, 0.0, 1.0)
    obs_sum = (obs * wgts).resample(time="AS").sum(dim="time")
    ones_out = (ones * wgts).resample(time="AS").sum(dim="time")
    obs_s = (obs_sum / ones_out).mean(("lat", "lon")).to_series()
    return obs_s

obs_s = weighted_temporal_mean(ds)

Convert the dataset into a dataframe

obs_df = pd.DataFrame(obs_s).rename(columns={0:'value'})

obs_df
value
time
1850-01-01 -0.338822
1851-01-01 -0.245482
1852-01-01 -0.291014
1853-01-01 -0.342457
1854-01-01 -0.276820
... ...
2017-01-01 0.777498
2018-01-01 0.641953
2019-01-01 0.809306
2020-01-01 0.865248
2021-01-01 0.683412

172 rows × 1 columns

Plot the Output

In this this case, we use hvPlot to plot the output!

smbb_plot = (anomaly_smbb.hvplot.scatter('time', height=300, width=500, color='lightgrey', legend=False, xlabel='Year', ylabel='Global Mean \n Temperature Anomaly (K)', hover=False) *\
anomaly_smbb.mean(1).hvplot.line('time', color='k', line_width=3, label='ensemble mean')  *\
obs_df.hvplot.line('time', color='red', label='observations')).opts(title='Smoothed Biomass Burning')


cmip6_plot = anomaly_cmip6.hvplot.scatter('time', height=300, width=500, color='lightgrey', legend=False, xlabel='Year', ylabel='Global Mean \n Temperature Anomaly (K)', hover=False) *\
anomaly_cmip6.mean(1).hvplot.line('time', color='k', line_width=3, label='ensemble mean') *\
obs_df.hvplot.line('time', color='red', label='observations').opts(title='CMIP6 Biomass Burning')

smbb_plot + cmip6_plot

Spin Down the Cluster

After we are done, we can spin down our cluster

cluster.close()
client.close()